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1.  Introduction 

Consider  tbs  system  of  llnssr  equations 
(1.1)  A x » b 


where  A Is  in  N«N  sparse  syznetrlc  positive  definite  mstrlx  such  as  those  that  arise  In 

finite  difference  and  finite  element  approximations  to  elliptic  boundary  value  problems 

In  two  and  three  dimensions.  A classic  method  for  solving  such  systems  Is  Gsusslan 

elimination:  We  use  the  k th  equation  to  eliminate  the  k th  variable  from  the  remaining 

N-k  equations  for  k • 1,2,  . . .,N-1  and  then  back-aolve  the  resulting  upper  triangular 

T 

system  for  the  unknown  vector  x.  Equivalently,  we  form  the  U DU  decomposition  of  A and 
successively  solve  the  triangular  systems 


(1.2) 


Uz-b,  Dy-z,  Ux-y 


Unfortunately,  as  the  elimination  proceeds,  coefficients  thst  were  zero  In  the  original 
system  of  equations  become  nonzero  (or  flll-ln) , Increasing  the  work  and  storage  required. 
The  purpose  of  this  r *P«r  1*  to  Introduce  a new  graph-theoretic  model  of  such  flll-ln; 
to  establish  lower  bounds  for  the  work  and  storage  associated  with  Gaussian  elimination; 
and  to  present  a minimal  storage  sparse  elimination  algorithm  that  significantly  reduces 
the  storage  required. 

In  section  2,  we  review  the  graph-theoretic  elimination  model  of  Farter  and  Rose 

and  Introduce  a new  element  model  of  the  elimination  process.  In  section  3,  we  use  this 

model  to  give  simple  proofs  of  Inherent  lower  bounds  for  the  work  end  storage  associated 

with  Gaussian  elimination,  generalizing  similar  reoults  of  Georgs  and  Hoffman,  Martin, 

and  Rose.  Last,  In  section  4,  we  show  how  the  element  model,  combined  with  the  rather 

unusual  Idea  of  recomputing  rather  than  saving  tha  factorization,  leads  to  a minimal 

storage  sparse  elimination  algorithm  that  requires  significantly  less  storage  than  regular 

2 2 

sparse  elimination  (e.g.  0(n  ) vs.  0(n  log  n)  for  the  five-  or  nine-point  operator  on  an 
n*n  mesh). 
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2.  An  Element  Model  for  Gaueeian  Elimination 


In  this  section,  we  shell  revlewr the  graph-theoretic  elimination  model  of  Gaussian 


elimination  "suggested  by  Parter  [6]  and  extensively  developed  by  Rose  [7]  and**lntroduce 
a new  element  model*,  cf.  Elsenstat  [1],  Sherman  [9]. 

Given  an  Irreducible  N»N  symmetric  positive  definite  matrix  A ■ (a^),  we  can 
represent  the  zero-nonzero  structure  of  A by  a graph  G(A)  ■ (V,E)  as  follows:  The  vertex 

v^  of  G(A)  corresponds  to  the  1 th  row/column  of  A;  the  edge  (v^.v^)  Is  an  edge  of  G(A)  If 
and  only  If  a^  i 0.  He  shall  model  Gaussian  elimination  on  the  matrix  A as  a sequence 
of  such  graphs. 

Let  » G(A).  Then  the  graph  is  derived  from  G^  as  follows:  Corre- 

sponding to  using  the  k th  equation  to  eliminate  the  k th  variable  from  the  remaining  N-k 
equations,  we  add  any  edges  necessary  to  make  all  vertices  adjacent  to  the  k th  vertex  v^ 
pairwise  adjacent  and  then  delete  v^  and  all  edges  Incident  to  It.  Thus  the  graph  G^ 
represents  the  nonzero  structure  of  the  lower  H-fc+1  * ff-k+1  submatrix  of  A Just  before  Che 
k th  variable  is  eliminated.  From  this  observation,  the  following  operation  and  storage 
counts  are  lamedlate: 


Theorem  (Rose  [7]):  Let  denote  the  degree  of  vertex  v^  In  the  elimination  graph  G^, 

l.e.  st  the  time  It  was  allmlnatad.  Then  the  number  of  multiplies  required  to  form  the 
UTDU  decomposition  of  A la  given  by 


* ■ kf1  f w» 


and  the  number  of  of f-dlegooal  nonzero  entries  In  U Is  given  by 


* - X 4.  . 

k-1 


The  element  model  emuletes  Geuaalan  elimination  as  e sequence  of  transformations 
on  the  greph  C(A)  and  the  collection  C of  maximal  cliques  of  vartlces  In  G(A)  which  we 


-'■rndtiv 


I 
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■hall  refer  to  u elmente.  Initially,  G*1*  - C(A) , E(1*  - E,  and  all  vartlcea  are 

narked  unellalnated.  Then  and  E^"*^  are  derived  from  G^  and  E^  aa  follows: 

Corresponding  to  ualng  the  kth  equation  to  ellnlnate  the  kth  variable  fron  the  renalnlng 

M-k  equations,  we  nark  th  kth  vertex  v^  aa  ellnlnated  and  add  any  edgea  necessary  to  sake 

all  vertices  adjacent  to  v^  pairwise  adjacent  (no  vertices  or  edges  are  deleted);  all  the 

elements  of  E^  containing  v^  are  merged  Into  a new  element  and  then  deleted  from  E^k\ 

Mote  that  the  subgraph  of  Induced  by  the  unellmlnated  vertices  la  Just  G^  so  chat 

the  number  of  unellmlnated  vertices  adjacent  to  v^  at  the  time  It  la  eliminated  Is  d^  as 

(k)  *(k) 

above.  Also,  elements  In  E are  cliques  in  C (though  not  necessarily  maximal 

-(k) 

cliques).  A vertex  will  be  said  to  be  an  exterior  vertex  In  C If  It  belongs  to  more 

than  one  clement;  otherwise  It  belongs  to  exactly  one  element  and  Is  said  to  be  an 
interior  vertex. 

A a an  example,  consider  the  elimination  process  for  the  nine-point  operator  on 
a 3*3  grid: 

1-5-2 

I I I 

Figure  2.1  7-8-9 

III 

3-8-4 

Initially,  the  elements  correspond  to  the  Individual  mesh  squares  or  elements,  since  these 
form  the  only  maxlsial  cliques  in  G(A).  Hence  the  nmse  "element  modal"  (cf.  George  [4]). 

When  we  eliminate  vertex  v^,  there  la  no  change  in  the  element  graph  (except  to  mark  v^ 
eliminated)  since  It  Is  an  Interior  vertex.  The  same  la  true  whan  we  eliminate  Vj,  Vj, 
and  v^.  When  we  eliminate  Vy  however,  the  two  elements  containing  Vj  are  merged,  as  also 
happens  when  we  eliminate  v^.  Mow  eliminating  v ^ causes  the  two  remaining  elements  to  be 
merged  Into  the  final  element  which  consists  of  all  the  vertices.  Eliminating  vertices 
Vg  and  v9  haa  no  further  affect  since  they  ara  Interior  vartlcea. 

Wa  note  the  following  properties  of  the  element  model  for  future  use.  Throughout 
the  elimination  procaaa,  the  number  of  elements  that  contain  a given  vertex  never  lncreaaes. 
Thus  an  eliminated  vertex  la  necessarily  an  Interior  vertex,  and,  correspondingly,  an 


4 

exterior  vertex  le  unellmlnated.  (However,  an  Interior  vertex  need  not  have  been 
eliminated.)  Initially,  the  largest  element  contalne  at  moat  d vertlcea,  where  d 

■IX  MX 

la  the  maximum  degree  of  any  vertex  In  G(A)  and  no  vertex  belonga  to  more  than  2<*max 

l(d_ax-l) 

element a.  (See  Eleenetat  [1]  for  an  example  where  a vertex  actually  belonga  to  22 
elemente.)  At  the  end  of  the  elimination,  the  only  element  conaleta  of  all  the  vertlcea. 

3.  Lever  Bounds  for  Gaussian  Elimination 

In  recent  yeare,  there  hea  been  a great  deal  of  Interest  In  computational  complexity, 
particularly  In  analysing  the  number  of  f undmaental  operations  Inherent  In  a computalon. 
Along  these  lines,  George  {4]  and  Hoffman,  Martin,  and  Rose  [5]  have  used  the  elimination 
graph  model  to  prove  lower  bounds  for  the  work  and  storage  associated  with  Gaussian 
elimination  on  symmetric  poaltlve  definite  matrlcae  whose  graphs  are  certain  regular 
planar  grids.  In  this  section,  we  ehall  describe  how  the  element  model  can  be  used  to 
give  simpler  proofs  of  these  results  that  extend  easily  to  some  Irregular  grids  In  two  and 
three  dimensions.  For  further  details,  sea  Elsenstat  (1]. 

First  we  shall  define  the  model  of  computation  In  which  these  lower  bounds  are 
valid.  Given  an  Irreducible  N«H  ay  metric  positive  definite  matrix  A,  we  seek  to  solve 
the  system  of  linear  equations  (1.1)  using  Gaussian  elimination.  The  work  associated  with 
the  elimination  process  will  be  taxon  as  the  number  of  multiplies  to  factor  A;  the  storage 
aa  the  ntmiber  of  nonsero  entries  In  the  matrix  U. 

The  work  and  storage  required  for  Gaussian  elimination  Is  directly  related  to  the 
order  In  which  the  variables  are  eliminated.  Thus , Instead  of  solving  the  original  system 
(1.1),  we  might  prefer  to  solve  the  permuted  system. 

3.1)  r A PT  y - P b,  PT  y - x 

for  soma  given  permutation  matrix  P so  aa  to  reduce  the  associated  work  and  storage. 

Indeed,  a great  deal  of  research  has  been  done  toward  discovering  good  ordering  strategies 
(e.g.  minimum  degree  [7]  and  neeted  dissection  [4]).  Mote  that  the  permutation  P does  not 


s 


change  tha  structure  of  tha  graph  C(A) , other  than  ralaballlng  the  vertices  and  thus 
changing  the  order  of  elimination.  Since  the  lower  bounds  we  shall  derive  will  be 
Independent  of  such  considerations,  they  will  be  valid  for  all  possible  orderings. 

As  In  section  2,  let  G(A)  ■ (V,E)  denote  the  graph  associated  with  A,  and  let  dmax 
denote  the  maximum  degree  of  any  vertex  In  G(A).  We  shall  now  make  the  following  additional 
assumption  about  the  matrix  A or,  equivalently,  the  graph  C(A): 

"iBOperitnetria"  Inequality  There  exists  constants  K > 0 and  0 < a,  0 £ 1 such  that, 
given  any  subset  S of  V with  |s|  S 8 N,  we  have 

l»S|  * K |S|“ 

where 

)S  ■ {v  < S :gw  t V-S  such  that  v and  w are  adjacent  In  G(A)} 
la  the  boundary  of  S. 

This  assiaiptlon  Is  closely  related  to  the  classical  lsoperlmetrlc  Inequality,  which 
relates  the  area  A and  perimeter  P of  a plane  figure 

2 lA 

and  the  laovolumetrlc  Inequality,  which  relates  the  volume  V and  surface  area  S of  a 
three-dimensional  region 

S3  k 36»V2. 

Indeed,  proofs  of  this  property  for  particular  grids  follow  the  classical  proofs.  As  an 
exMple,  for  the  five-  or  nlne-oolnt  operator  on  an  n*n  grid,  we  have 

|»S|  * 1S| 1/2  If  |S|  S ^ , 

while  for  the  seven-  or  twenty-seven-point  operator  on  an  n*nxn  grid  we  have 
|»S|  l |S|2/3  If  ISJ  i ^ , 


I 


cf.  Elsenstat  [1J. 


6 

Lama  (Elsenatat  ( 1 ] ) ; Assume  that  tha  graph  G(A)  aatlaflaa  an  laoparlaatrlc  Inequality. 
Then  thera  exists  a K(8N/2<1,“J,)a  clique  In  aoae  elimination  graph  G^k\ 

Proof'.  At  the  atart  of  the  elimination  proceaa,  the  maximum  size  of  any  element  la  dggg, 
the  maxlmlmum  degree  of  any  vertex  In  G(A).  At  the  end  of  the  elimination  proceaa,  there 
la  exactly  one  element  of  alze  N,  namely  the  act  of  all  vertlcea.  Therefore,  at  acme  point 
during  the  elimination,  the  flrat  element  of  size  > 6N  was  created.  This  element  was 
created  by  merging  all  the  elements  containing  some  vertex  v^.  As  we  saw  In  section  2, 
could  belong  to  at  most  l1*01**  elements  at  this  point  and  thus  at  least  one  such,  say  e*, 
must  contain  mote  than  BN/2<*,“X  vertices.  Consider  the  boundary  vertices  of  e*.  There 
are  at  least  KCSN/l**®**)0  such  vertices  by  the  isoperlmetrlc  Inequality,  all  are  uneltml- 
nated  vertices  at  this  stage,  and  they  are  all  pairwise  adjacent.  Thus  they  form  a clique 
In  C(k).  QED 

Corollary : The  bandwidth  of  the  matrix  A la  at  least 

K(B/2d»“)a  Na 

Corollary : The  total  work,  required  to  factor  the  matrix  A using  band  elimination  Is  at 

least 

~ *2(fl/2<lB"t)2ct  N(2a+1)  ; 

the  nuofcer  of  nonzeroes  In  the  band  of  A la  at  least 
C(e/2dmax)01  jj(«*1)  t 

The  proof  follows  directly  from  the  standard  operation  and  storage  counts  for  band 
elimination  In  terms  of  the  number  of  equations  and  the  bandwidth. 

Corollary:  The  total  work  required  to  factor  the  matrix  A is  at  least 

| K3(*/2d*»*)3a  N3®; 


the  number  of  nonxeroes  In  U la  at  least 
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i K2(6/2d®**)2a  N2°. 

The  proof  follows  fro*  the  Lesna  as  In  Hoffman,  Martin,  and  Rose  [5]. 

We  shall  now  examine  the  consequences  of  these  simple  results.  For  planar  grids 
2 

with  ~n  vertices  and  bounded  degree,  such  as  the  five-  or  nine-point  operator  on  an  n«n 
grid: 

(1)  the  bandwidth  Is  at  least  0(n) 

4 3 

(2)  the  work  and  storage  are  at  least  0(n  ) and  0(n  ) respectively  for  band  elimination 

3 2 

(3)  the  work  and  storage  are  at  least  0(n  ) and  0(n  ) respectively  for  sparse 

2 

elimination  (a  more  careful  analysis  gives  0(n  log  n)  for  the  storage;  see 
Elsenstat  (1]). 

For  three  dimensional  grids  with  ~n^  vertices  and  bounded  degree,  such  as  the  seven-  or 
twenty-seven-point  operator  on  an  n*n«n  grid: 

2 

(4)  the  bandwidth  is  at  least  0(n  ) 

(5)  the  work  and  storage  are  at  least  0(n2)  and  O(n^)  respectively  for  band  elimination 

6 4 

(6)  the  work  and  storage  are  at  least  0(a)  and  0(n  ) respectively  tor  sparse 
elimination. 

4.  Minimal  Storage  Sparge  Elimination 

One  of  the  major  disadvantages  of  Gaussian  elimination  for  solving  sparse  systems 

of  linear  equations  Is  the  amount  of  storage  required.  For  example,  to  solve  the  nine-point 

2 

finite  difference  operator  on  an  n*n  grid  requires  at  least  0(n  log  n)  storage  whereas  an 

2 

Iterative  method  would  require  only  0(n  ) storage.  In  this  section,  we  present  a variation 
of  sparse  Gaussian  elimination  that  tiades  a significant  reduction  In  storage  for  a modest 
Increase  In  work.  For  ease  of  exposition,  we  shall  restrict  attention  to  the  nine-point 
operator  on  an  n«n  grid  with  n * 2C-1,  but  the  results  are  valid  for  more  general  finite 
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element  grids,  cf.  Elsenstat,  Schultz,  and  She  naan  [3],  Sherman  [9].  For  similar  results 
on  band  elimination,  see  Elsenstat,  Schultz,  and  Sherman  [2],  Sherman  [9]. 

Two  baalc  concepts  are  uaed  In  achieving  this  reduction:  First,  rather  than  save 
the  entire  factorization,  we  shall  threw  most  of  the  nonzero  entries  of  U away  and  re- 
compute them  as  needed  during  the  back-solution;  second,  ve  use  an  element  merge  tree  to 
specify  a dlvlde-and-conquer  elimination  ordering  and  to  keep  track  of  those  entries  of  U 
that  are  being  saved  during  each  step  of  the  calculation. 

Suppose  we  were  to  perform  Gaussian  elimination  In  such  a way  that  the  last  2n+l 
variables  to  be  eliminated  correspond  to  the  vertices  on  a dividing  cross  as  shown  In 
Figure  4.1.  At  the  end  of  the  elimination  process,  we  would  have  generated  all  the 
coefficients  In  the  upper  triangular  system  of  equations  that  remains,  and  we  have  merely 
to  solve  this  system  for  the  vector  of  unknowns  x.  Yet  suppose  we  had  saved  only  those 
coefficients  in  the  last  2n+l  rows.  Then  we  could  only  solve  for  the  last  2n+l  variables, 
l.e.  the  values  of  the  unknowns  on  the  dividing  cross.  This  Is  enough,  however,  to  split 
our  original  n*n  problem  into  four  smaller  ~n/2  * ~n/2  problems  of  the  same  form,  which  we 
can  now  solve  in  the  same  fashion. 

I 

Figure  4.1. 

Ill 

Of  course,  we  will  have  to  do  more  work  than  we  would  have  had  we  saved  the  entire 
factorization.  But  how  much  more?  The  nested  dissection  ordering  of  George  [4]  orders 
the  variables  on  the  dividing  cross  last  and  requires  approximately  Cn^  multiplies  for 
gome  fixed  constant  C.  Thus  the  cost  of  the  whole  procedure  Is  Just  Cn^  plus  the  cost  of 
solving  four  ~n/2  * — n/2  problems.  Letting  8(n)  denote  this  cost,  we  have  that 

8(n)  » Cn^  ♦ 48(n/2) 

and  that  8(n)  - Cn"5  for  n sufficiently  small.  (After  all,  when  the  amount  of  storage  for 


the  factorization  gets  small  enough,  we  may  as  well  save  all  of  It.)  The  solution  to  this 
recurrence  relation  Is 

6(n)  » 2Cn2 3. 

Thus  we  are  doing  twice  as  much  work  but,  as  we  shall  see,  the  savings  In  storage  will  be 
much  more  significant. 

The  element  merge  tree  is  based  on  the  element  model  Introduced  in  section  2. 

Recall  that  Initially  there  were  a number  of  elements,  or  maximal  cliques  of  vertices,  and 
Chat  as  the  elimination  progressed  these  elements  were  merged  into  larger  and  larger 
elements  until  only  a single  element  containing  Che  entire  set  of  vertices  remained. 

Given  an  elimination  ordering,  we  construct  the  corresponding  element  merge  tree  as  follows 
The  nodes  In  the  tree  represent  elements  that  were  created  during  the  elimination  process; 
the  root  of  the  tree  Is  the  final  element  consisting  of  the  entire  set  of  vertices;  the 
node  (l.e.  element)  e^  Is  a son  of  another  node  e^  If  and  only  if  e^  was  merged  Into  e^ 
when  some  vertex  was  eliminated.  The  merge  tree  for  the  3*3  nine-point  operator  of 
Figure  4.2a  Is  given  In  Figure  4.2b. 


1-5-2 

I I i 

Figure  4.2a.  7-8-9 

I I I 

3-6-4 


Figure  4.2b. 


(1,2,3,. ...9) 


(1,2, 5. 7, 8, 9) 

, / \ 

(1, 5,7,8)  (2,5, 8, 9) 


(3, 4, 6, 7, 8, 9) 

/ \ 

(3,6, 7.8)  (4,6, 8,9) 


We  will  now  specify  the  ordering  of  vertices  and  the  corresponding  element  tree 

2 

which  we  trill  use  to  make  the  storage  required  for  our  procedure  0(n  ) as  opposed  to  the 

2 

0(n  log  n)  for  the  factorization.  Beginning  with  the  entire  grid,  we  break  It  into  four 
equal-size  elements  using  a dividing  cross  as  In  Figure  4.1.  (Note  that  the  vertices  on 
the  dividing  cross  belong  to  more  than  one  element  since  they  are  exterior  vertices.)  The 
last  vertices  to  be  eliminated  will  be  the  center  vertex  followed  by  the  other  vertices  on 


the  dividing  cross.  Ue  then  order  the  Interior  vertices  In  element  I,  followed  by  those 
of  element  II,  element  III,  and  element  IV.  The  vertices  within  each  of  these  are 
eliminated  In  an  analogous  fashion:  The  last  vertices  to  be  eliminated  are  the  center 

vertex  followed  by  the  other  vertices  on  a dividing  cross;  we  then  order  the  vertices  In 
each  of  the  four  subelements,  and  so  on. 

The  ordering  that  results  can  be  shown  to  be  equivalent  to  the  nested  dissection 
ordering  In  terms  of  the  work  and  storage  required  (cf.  Elsenstat,  Schultz,  and  Sherman  [3], 
Rose  and  Whitten  [8],  Sherman  [9]),  and  the  element  merge  trees  are  Identical,  although  the 
actual  orderings  are  completely  different.  Thus  the  divide- and  conquer  order  produces  an 

O(n^)  elimination  scheme.  We  now  show  that  the  elimination  can  be  carried  out  In  such  a 

2 

way  thst  the  coefficients  in  the  last  2n+l  rows  can  be  computed  using  only  0(n  ) storage. 
Consider  again  the  element  merge  tree.  The  nodes  at  the  bottom  or  zeroth  level 

0 2 

correspond  to  elements  in  the  original  graph  G(A),  all  of  which  contain  exactly  4 « (2  +1) 
vertices.  At  the  first  level,  the  elements  were  formed  by  merging  four  bottom-level  elements 
and  thus  each  contains  exactly  9 ■ (2  +1)  vertices.  In  general,  at  the  k th  level  the 

elements  were  formed  by  merging  four  elements  on  the  (k-l)t?i  level  and,  by  induction, 
k 2 

contain  precisely  (2  +1)  vertices. 

We  shall  say  that  a particular  element  in  the  element  tree  is  active  at  a particular 
point  in  the  elimination  if  It  has  not  yet  been  merged  into  another  element.  Then  the  only 
entries  in  the  triangular  factor  U that  we  shall  be  saving  at  that  point  are  those  corre- 
sponding to  pairs  of  exterior  vertices  in  active  elements.  Note  that  the  final  element 
Is  active  once  It  Is  created  and  Is  never  deactivated,  so  that  the  coefficients  In  the  last 
2n+l  rows  of  U will  be  available  to  solve  for  the  unknowns  on  the  dividing  cross  at  the 
end  of  the  elimination. 

Since  there  are  at  most  four  active  elements  at  any  level  in  the  element  merge  tree 
(other  than  the  zeroth)  at  any  stage  of  the  elimination  by  virtue  of  the  recursive  definition 
of  the  ordering,  we  cen  easily  bound  the  total  storage  required:  At  the  k th  level  from  the 

bottom,  there  are  at  most  1 ■ 4"2k  exterior  vertices  per  element,  which  require  at  most 
1(1+1) /2  nonzero  entries  of  U to  be  saved;  thus  the  total  storage  is  at  most 
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V 4 • 4 (*-2k)  (4-2k+l)  » n2. 

k-0  2 

A tore  careful  analysis  using  two-way  dissection  show  that  the  total  storage  can  be  reduced 

to  9n2,  cf.  Elsenstat , Schultz,  and  Sherman  {3],  Sherman  [9]. 

' Note  that  this  surprising  result  does  not  really  violate  the  results  of  section  3; 

2 

the  storage  for  U still  la  0(n  log  n)  — we  Just  aren't  saving  all  the  entries.  Another 
point  worth  noting  Is  that  the  algorithm  can  be  Implemented  to  run  In  time  0(n2°82  2)  If 
the  elimination  of  variables  on  dividing  crosses  Is  done  In  blocks  using  Strassen's 
algorithm  for  matrix  multiplication  and  Inversion  [10].  Again,  this  does  not  violate 
the  results  of  section  3 since  we  are  doing  block  rather  than  point  elimination. 
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